CHC analysis
Dan Bruzzese
6/17/2020
using chc 6
present CV LDAs, have testing/training if need be
graphed LDAs are all the points.
MANOVA for species level
I cannot do site contrasts because of our uneven samplign scheme. The MANOVA throws out samples.
fit_full_species_man$pca_summary ## Importance of first k=7 (out of 35) components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.0736 0.7112 0.47335 0.37674 0.33980 0.30966 0.2955
## Proportion of Variance 0.4066 0.1784 0.07903 0.05006 0.04073 0.03382 0.0308
## Cumulative Proportion 0.4066 0.5850 0.66402 0.71409 0.75481 0.78864 0.8194summary(manova(as.matrix(data[,4:cols]) ~ hostRace * sex, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## hostRace 4 0.03980 40.088 28 776.62 < 2.2e-16 ***
## sex 1 0.81535 6.956 7 215.00 1.795e-07 ***
## hostRace:sex 4 0.79391 1.833 28 776.62 0.005633 **
## Residuals 221
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PCA for Species Level
LDA for Species level
## Reference
## Prediction Cingulata Cornivora Mendax pom Zepheria
## Cingulata 25 0 0 0 0
## Cornivora 1 5 0 0 0
## Mendax 0 0 60 5 1
## pom 0 0 1 116 0
## Zepheria 0 0 0 1 16## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 9.610390e-01 9.390180e-01 9.273307e-01 9.820323e-01 5.281385e-01
## AccuracyPValue McnemarPValue
## 1.508463e-49 NaN
Bayesian Model
Collapsed all 35 chc measurements into a value column with a corresponding chc column detailing what chc was measured. Fitted a model to evaluate the efect of host and sex on chc values. The model accounts for uneven site effects and individual effects.
All Bayesian models were created in Stan computational framework (http://mc-stan.org/) accessed with brms package (Bürkner, 2017). To improve convergence and guard against overfitting, we specified mildly informative conservative priors. We ran 4 chains with 2000 interations. We checked for convergence and found all models mixed well. (Pareto K estimates are all good!)
#1 value ~ hostRace * sex + (1|site) +(1|chc) + (1|ID)
model comparison?
Model Converged
mod = readRDS("fit_chc_model1.rds")
plot(mod)
model fit is good = posterior values match real values
brms::pp_check(mod, resp ="value") #great## Using 10 posterior samples for ppc type 'dens_overlay' by default.
Model fit is good = r2 == 88%
bayes_R2(mod, summary = TRUE)## Estimate Est.Error Q2.5 Q97.5
## R2 0.8790622 0.0009515364 0.8771306 0.8808946model summary
tab_model(mod)| Â | value | |
|---|---|---|
| Predictors | Estimates | CI (95%) |
| Intercept | -0.25 | -0.49 – -0.00 |
| hostRace: Cornivora | -0.13 | -0.29 – 0.03 |
| hostRace: Mendax | 0.17 | 0.09 – 0.26 |
| hostRace: Zepheria | 0.33 | 0.23 – 0.43 |
| hostRace: pom | 0.11 | 0.04 – 0.17 |
| sex: Male | 0.03 | -0.06 – 0.11 |
| hostRaceCornivora.sexMale | -0.14 | -0.35 – 0.06 |
| hostRaceMendax.sexMale | -0.07 | -0.17 – 0.03 |
| hostRaceZepheria.sexMale | -0.01 | -0.14 – 0.12 |
| hostRacepom.sexMale | -0.14 | -0.23 – -0.05 |
| Random Effects | ||
| σ2 | 0.44 | |
| Ï„00 | 0.07 | |
| ICC | 0.86 | |
| N site | 8 | |
| N ID | 231 | |
| N chc | 35 | |
| Observations | 8085 | |
| Marginal R2 / Conditional R2 | 0.021 / 0.879 | |
two way interactions
plot_model(mod, type = "pred", terms = c("hostRace", "sex"))
fixed effects
plot_model(mod, bpe = "mean", bpe.style = "dot", prob.inner = 0.4, prob.outer = 0.8)
hypothesis testing?
Manova for pomonella
summary(manova(as.matrix(data[,4:cols]) ~ hostRace * site * sex, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## hostRace 1 0.82295 3.3194 7 108 0.003071 **
## site 1 0.85799 2.5537 7 108 0.017937 *
## sex 1 0.52503 13.9576 7 108 8.859e-13 ***
## hostRace:site 1 0.69396 6.8041 7 108 1.085e-06 ***
## hostRace:sex 1 0.95188 0.7799 7 108 0.605479
## site:sex 1 0.93727 1.0326 7 108 0.412735
## hostRace:site:sex 1 0.90238 1.6691 7 108 0.124190
## Residuals 114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1LDA for Pomonella
## Reference
## Prediction Apple Haw
## Apple 25 12
## Haw 14 71## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.786885246 0.503288443 0.703538181 0.855813669 0.680327869
## AccuracyPValue McnemarPValue
## 0.006229331 0.844519267
LDA predicting host and sex for pomonella
## Reference
## Prediction Apple_Female Apple_Male Haw_Female Haw_Male
## Apple_Female 23 2 9 1
## Apple_Male 1 1 1 2
## Haw_Female 8 0 32 3
## Haw_Male 1 3 5 30## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 7.049180e-01 5.725131e-01 6.155825e-01 7.840219e-01 3.852459e-01
## AccuracyPValue McnemarPValue
## 8.606186e-13 9.110307e-01
MANOVA for Zepheria
# no interaction because missing males at MtPleasant
summary(manova(as.matrix(data[,4:cols]) ~ sex + site, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.74708 0.67708 5 10 0.6507
## site 1 0.83749 0.38808 5 10 0.8461
## Residuals 14LDA for Zepheria
## Reference
## Prediction Zepheria_Female Zepheria_Male
## Zepheria_Female 5 3
## Zepheria_Male 4 5## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.5882353 0.1793103 0.3292472 0.8155630 0.5294118
## AccuracyPValue McnemarPValue
## 0.4062810 1.0000000
LDA for sex +site zepheria
## Reference
## Prediction Zepheria_Female_EastLansing
## Zepheria_Female_EastLansing 2
## Zepheria_Female_MtPleasant 0
## Zepheria_Male_EastLansing 4
## Reference
## Prediction Zepheria_Female_MtPleasant
## Zepheria_Female_EastLansing 2
## Zepheria_Female_MtPleasant 1
## Zepheria_Male_EastLansing 0
## Reference
## Prediction Zepheria_Male_EastLansing
## Zepheria_Female_EastLansing 1
## Zepheria_Female_MtPleasant 5
## Zepheria_Male_EastLansing 2## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.29411765 -0.05699482 0.10313551 0.55958272 0.47058824
## AccuracyPValue McnemarPValue
## 0.95792652 0.03207164
MANOVA for Mendax
summary(manova(as.matrix(data[,4:cols]) ~ sex * site, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.57727 6.1024 6 50 7.541e-05 ***
## site 2 0.10716 17.1232 12 100 < 2.2e-16 ***
## sex:site 2 0.74705 1.3082 12 100 0.2258
## Residuals 55
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1LDA for Mendax
## Reference
## Prediction Mendax_Female Mendax_Male
## Mendax_Female 20 15
## Mendax_Male 13 13## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.54098361 0.07072905 0.40849889 0.66935590 0.54098361
## AccuracyPValue McnemarPValue
## 0.55241652 0.85010674
LDA for sex +site mendax
## Reference
## Prediction Mendax_Female_Fenville Mendax_Female_OtisLake
## Mendax_Female_Fenville 7 1
## Mendax_Female_OtisLake 1 5
## Mendax_Female_Sewanee 0 2
## Mendax_Male_Fenville 6 0
## Mendax_Male_OtisLake 0 1
## Mendax_Male_Sewanee 0 0
## Reference
## Prediction Mendax_Female_Sewanee Mendax_Male_Fenville
## Mendax_Female_Fenville 0 5
## Mendax_Female_OtisLake 1 0
## Mendax_Female_Sewanee 4 1
## Mendax_Male_Fenville 0 5
## Mendax_Male_OtisLake 1 1
## Mendax_Male_Sewanee 4 0
## Reference
## Prediction Mendax_Male_OtisLake Mendax_Male_Sewanee
## Mendax_Female_Fenville 0 0
## Mendax_Female_OtisLake 1 2
## Mendax_Female_Sewanee 2 3
## Mendax_Male_Fenville 0 1
## Mendax_Male_OtisLake 1 1
## Mendax_Male_Sewanee 2 3## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.409836066 0.283523654 0.285504382 0.543223627 0.229508197
## AccuracyPValue McnemarPValue
## 0.001280557 NaN
MANOVA for cingulata
summary(manova(as.matrix(data[,4:cols]) ~ sex * site, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.78806 1.2774 4 19 0.3136059
## site 1 0.38209 7.6817 4 19 0.0007371 ***
## sex:site 1 0.70133 2.0229 4 19 0.1319094
## Residuals 22
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1LDA for Cingulata
## Reference
## Prediction Cingulata_Female Cingulata_Male
## Cingulata_Female 5 7
## Cingulata_Male 7 7## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.46153846 -0.08333333 0.26587122 0.66629178 0.53846154
## AccuracyPValue McnemarPValue
## 0.83731354 1.00000000
LDA for sex +site cingulata
## Reference
## Prediction Cingulata_Female_SouthBend Cingulata_Female_Urbana
## Cingulata_Female_SouthBend 1 2
## Cingulata_Female_Urbana 2 0
## Cingulata_Male_SouthBend 4 0
## Cingulata_Male_Urbana 3 0
## Reference
## Prediction Cingulata_Male_SouthBend Cingulata_Male_Urbana
## Cingulata_Female_SouthBend 5 2
## Cingulata_Female_Urbana 4 1
## Cingulata_Male_SouthBend 0 0
## Cingulata_Male_Urbana 1 1## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.076923077 -0.243027888 0.009455391 0.251302917 0.384615385
## AccuracyPValue McnemarPValue
## 0.999943167 0.389256527
MANOVA for Cornivora
There is only one site
# no site because only sampled at one site.
summary(manova(as.matrix(data[,4:cols]) ~ sex, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.54542 0.83346 2 2 0.5454
## Residuals 3LDA for cornivora
Cant show a plot
## Reference
## Prediction Cornivora_Female Cornivora_Male
## Cornivora_Female 2 1
## Cornivora_Male 0 2## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.8000000 0.6153846 0.2835821 0.9949492 0.6000000
## AccuracyPValue McnemarPValue
## 0.3369600 1.0000000